Genetic tracing of the illegal trade of the white-bellied pangolin (Phataginus tricuspis) in western Central Africa

The white-bellied pangolin is subject to intense trafficking, feeding both local and international trade networks. In order to assess its population genetics and trace its domestic trade, we genotyped 562 pangolins from local to large bushmeat markets in western central Africa. We show that the two lineages described from the study region (WCA and Gab) were overlapping in ranges, with limited introgression in southern Cameroon. There was a lack of genetic differentiation across WCA and a significant signature of isolation-by-distance possibly due to unsuspected dispersal capacities involving a Wahlund effect. We detected a c. 74.1–82.5% decline in the effective population size of WCA during the Middle Holocene. Private allele frequency tracing approach indicated up to 600 km sourcing distance by large urban markets from Cameroon, including Equatorial Guinea. The 20 species-specific microsatellite loci provided individual-level genotyping resolution and should be considered as valuable resources for future forensic applications. Because admixture was detected between lineages, we recommend a multi-locus approach for tracing the pangolin trade. The Yaoundé market was the main hub of the trade in the region, and thus should receive specific monitoring to mitigate pangolins’ domestic trafficking. Our study also highlighted the weak implementation of CITES regulations at European borders.


Mitochondrial DNA
The NJ phylogenetic tree based on 641 cyt b sequences supported the structuring into six WBP geographic lineages, including Western Africa, Ghana, Dahomey Gap, WCA, Gab and Central Africa (Supplementary Fig. S2 online).A total of 579 individuals from Cameroon, Equatorial Guinea and northern Gabon (Bitam) clustered within WCA, whereas ten samples from southern and northern Gabon and three samples from central-south Cameroon (Yaoundé and Abong Mbang) clustered within Gab (Fig. 1).
A total of 67 and eight haplotypes were identified in WCA and Gab, respectively.WCA and Gab were differentiated by 23 mutations in the median-joining haplotype network, whereas no particular geographic structuring could be observed within the two lineages (Supplementary Fig. S4 online).Three haplotypes were dominant in WCA, with H2 (28%) and H3 (25%) located only in Cameroon, and H5 (17%) in Cameroon and Bioko Island (Supplementary Fig. S5 online).In Gab, the dominant haplotype (H70; 50%) was distributed from northern to southeastern Gabon.Douala markets recruited 8-17 haplotypes (dominated by H5 and H2 at 27% and 23%, respectively), while the Yaoundé market was characterized by 35 haplotypes (dominated by H3 at 32%).Haplotype (Hd) and nucleotide (π) diversity were respectively 0.830 and 0.010 in WCA, and 0.897 and 0.019 in Gab (Supplementary Table S2

Genetic structure
The PCoA did not show any clear structuring of the genetic space within WBP from WCA and Gab (Supplementary Figure S7 online).
Clustering analysis with STRU CTU RE, within all individuals and with prior information on lineage detected K = 2 as the most likely number of clusters (Fig. 3A).The two clusters did not totally correspond to the mitochondrial geographic assignment of the individuals regarding WCA and Gab lineages (Fig. 3A).In this case, 98% of all WBP individuals of WCA clustered into one population.Five out of eight mtDNA-assigned to Gab individuals clustered into a single group (Gab lineage).Three individuals (Y110, Ab10 and Ab2) from Yaoundé and Abong Mbang in Cameroon that were mtDNA-assigned to Gab clustered with WCA.Three individuals from Yaoundé and Douala (Y360, Y101 and DlaB103) that were mtDNA-assigned to WCA clustered with Gab.One individual from Yokadouma (Ou4) in Cameroon showed an admixed pattern between WCA and Gab (c.60/40% assignment, respectively).From K = 3, no specific population structure was found (Supplementary Fig. S8 online).
Within WCA and among reference populations, STRU CTU RE detected K = 2 as the most likely number of clusters (Fig. 3B).Two 'pure' genetic clusters regrouped Bayomen, Yabassi, and Foumbot from western and central Cameroon on one side, and PNCM, Abong Mbang and Sangmelima from North of Equatorial Guinea and southern / central Cameroon on the other side.Several reference populations from western, central and southern Cameroon, North of Gabon and South of Equatorial Guinea were admixed (Eseka, GES, Manengole and Mt Cam), and did not correspond to coherent geographic delineations (Fig. 3B).At K = 3, the four individuals from Medouneu (northern Gabon; GES population) clustered into an exclusive group.Seven individuals www.nature.com/scientificreports/from Sangmelima showed full to admixed assignment to the Medouneu cluster.From K = 4, no further coherent geographic structuring was found (Supplementary Fig. S9 online).Geneland analysis within WCA did not allow to geographically delineate consistent populations.Optimal number of populations varied between 10 and 13 clusters (Supplementary Fig. S10 online), with 11 as the most probable number of clusters according to the second-step iterative run (Supplementary Fig. S10 online).Assignment of individuals to populations among the best five runs was greatly unstable (data not shown).
We detected a significant IBD pattern among WBP individuals across the studied range (r = 0.099, p = 0.001), whereas there was no significant effect among populations (r = 0.298, p = 0.094; Supplementary Fig. S11 online).

Demographic history
We detected a significant signal of bottleneck in WCA with the S.M.M. model (Wilcoxon sign-rank test; p < 0.05), whereas the T.P.M. model did not (p = 0.061).
Tracing the pangolin trade A total of 535 (96%) of 557 WBP samples had unique genotypes.Eleven pairs of samples from Djoum, Manengole, Yaoundé and Medouneu shared a same genotype.The null hypothesis of encountering the same genotypes more than once by chance was rejected for all pairs of individuals (P < 0.0001).Unbiased probability of identity (uPI) and probability of identity among siblings (PIsibs) were both very low (uPI = 0.00125 ± 0.00538; PIsibs = 0.0236 ± 0.07236; Supplementary Fig. S12 online).
Seven of a total of 11 loci that presented appropriate private allele signatures for one or several populations using private allele rarefaction, could be related to 18 private alleles observed in the reference populations (WCA; Supplementary Table S6 online).On this basis, we could trace a total of 47 (12.8%) individuals sold as part of the www.nature.com/scientificreports/bushmeat trade in the study region (Supplementary Table S7 online) to seven potential reference populations (see Fig. 5).Thirty-five WBP from the Yaoundé market were harvested from Sangmelima, Eseka, Mt Cameroon, Abong Mbang, Manengole and South of Equatorial Guinea (GES).Ten individuals sold in Douala markets came from Abong Mbang, Sangmelima and Yabassi.Two pangolins seized at Brussels airport were traced to South of Equatorial Guinea and Abong Mbang.

White-bellied pangolins from western central Africa: range-overlapping lineages and weak geographic structuring
Our exhaustive sample set across three bordering countries, representing 563 newly collected samples of WBP from Cameroon, Equatorial Guinea and Gabon, allowed refining the lineage distribution of WBP from western central Africa.We show that the two mitochondrial lineages described from the study region are overlapping in ranges, with Gab being present all over Gabon and reaching southern Cameroon, while the southern extension of WCA reaches northern Gabon.Thus, our study clarifies the geographic distribution of the two lineages in western central Africa as previously based on a reduced number of samples 20 , and suggests that Gab is not endemic to Gabon.The three individuals representing Gab in Cameroon were found in the Yaoundé market and in the reference population of Abong Mbang.The latter corresponds to a local market that is off the main road network connecting northern Gabon to southern Cameroon 27 .The involved market stakeholder indicated that the two pangolins originated from the Lomié forest (c. 100 km from the selling site), thus suggesting that the Gab lineage naturally occurs in Cameroon.
Microsatellites genotyping showed that WCA and Gab corresponded to differentiated evolutionary units (mean F ST = 0.192), in line with previously published nDNA sequencing data 20 and with the mtDNA pattern found in the present study.However, because the three mtDNA-assigned Gab individuals from Cameroon were cyto-nuclear hybrids (i.e., assigned to the WCA cluster using microsatellites), ancient gene flow between the two lineages might be involved and is blurring their geographic delineation, notably in their supposed range of co-occurrence.All the more since an admixed individual from Yokadouma, in south-eastern Cameroon (Yokadouma), provides evidence for ongoing (but probably limited) introgression between the two lineages.Further sampling in forests from southern Cameroon and northern Gabon will have to be conducted in order to achieve a more comprehensive picture of the respective ranges of the WCA and Gab genetic entities, and their potential zones of admixture.
Despite the use of hypervariable microsatellites markers and a fair level of genetic diversity observed across populations ( 21 ; this study) along the c. 1,000 km extension of our study zone, genetic structuring and differentiation within WCA was weak.Structure was the only analysis capable of detecting population clusters, using the available signal of geographic structuring based on prior delineation of populations (r < 1; data not shown).However, the two 'pure' genetic clusters did not correspond to anything coherent in space, with a first cluster located strictly south of the Sanaga river (PNCM + Abong Mbang + Sangmelima) and the other spread on both sides of the river (Bayomen [north] vs. Yabassi and Foumbot [south]), while the admixed individuals were also present north and south of the river (Manengole and Mt Cam [north] vs. Eseka and GES [south]).www.nature.com/scientificreports/At K = 3, WBP from northern Gabon (Medouneu), which corresponds to the southernmost limit of WCA, were identified as a distinct population cluster.Because Medouneu is located c. 50 km from Monts de Cristal, a putative Pleistocene rainforest refugium 28 , pangolins from this area may hold the signature of a relict population, as previously observed in mammals and plants.However, the lack of population structuring in other areas of putative refugia such as southern Gabon and western Cameroon 28 , both sampled in our study, questions this interpretation; especially since Medouneu is located at the bottom end of the North-South IBD pattern that we observed in WCA (see below), which could be the cause of such apparent clustering (see 29 ).The refugium hypothesis is also challenged by the presence of Medouneu-related individuals (pure or admixed) further north in southern Cameroon (Sangmelima), suggesting a more extended range for this population.
The lack of genetic structure and population differentiation across WCA is to be interpreted in the light of the IBD pattern that we detected among WBP individuals.Although our study did not quantify dispersal rate and range, the pattern that we observed fits a stepwise migration model where dispersal distances are large enough to lead to near-panmictic genetic patterns 30 .Our results confirm the non-influential role of the Sanaga river on the biogeography and population delimitation of several African taxa (see 31 ), although its role as a limiting -i.e., not prohibitive-factor to gene flow (e.g., 32 ) was not tested here.They also suggest that there is no specific barrier to gene flow for WBP in the study zone (WCA), in line with the lack of structure found in another lineage of WBP investigated with the same markers (Dahomey Gap; 22 ).Our population genetic study again provides evidence for unsuspected dispersal capacities in WBP, which is also coherent with a recent, spatial and/or demographic expansion across the WCA range following forest post-glacial expansion during the Late Pleistocene (e.g., 33 ), as suggested by the mtDNA analyses.Despite habitat fragmentation and hunting activities occurring in the study zone 14 , dispersal may still be maintained given the capacity of the species to cross and occur in anthropized habitats 34 and the still prevalent, large forested areas available 35 .Field surveys are blatantly needed to document the dispersal behavior of WBP, notably during their early phase of supposed vagrancy 36 .
The only individual from Bioko Isl. that we could sequence and genotype did not show any genetic differentiation from the continent.Thus, the question of whether WBP sold on Bioko Isl. are imported from the continent (Cameroon) or not cannot be answered here 8 , and as a corollary, that of a surviving population endemic to the island.Additional genetic surveys in Bioko Isl., both from bushmeat markets and field sites, will have to be www.nature.com/scientificreports/conducted in order to assess the status of WBP on the island, an urgent matter of conservation as island populations are often subject to rapid genetic drift, accumulation of deleterious variants and inbreeding.

Highest levels of genetic diversity and Wahlund effect in the context of a recent demographic decline
Overall, populations from the WCA and Gab lineages were characterized by the highest levels of genetic diversity described in WBP.This is important for the conservation of the species in western central Africa, as fair levels of genetic diversity are expected to promote population fitness and ensure adaptability to environmental change 37 .MtDNA nucleotide diversity was by far the highest for WCA (largest sampling effort; N > 550) and Gab (reduced sampling effort, similar to other lineages; N = 13), c. 1.7-9.5 times greater compared to the other lineages.Mean allelic richness in WCA was also higher than in the Dahomey Gap lineage (A R = 10.79 vs. 4.27; minimum sample size = 169; 21,22 ), the only comparable study available.
The generally significant deficit of heterozygotes observed among populations (F IS = 0.112-0.217)could lead to the conclusion that WBP populations from WCA are inbred.However, inbreeding would come in contradiction with the fairly high levels of heterozygosity observed among populations (mean Ho = 0.579; mean He = 0.664) and the weak geographic structuring observed across WCA (see 38 ).Because population differentiation (F ST ) across WCA was low whereas F IS values were high, our results fit the expectation of a Wahlund effect, i.e.where a same population regroups individuals from genetically differentiated entities 39 .Such signature could be reached through regular, long-range inter-connectivity among WBP populations, implying again unsuspected dispersal capacities of the species.However, given the knowledge gaps affecting the natural history of WBP and African pangolins in general 4 , we cannot discard the fact that some unknown breeding strategy of the species could participate to this Wahlund effect.
We identified an important decline in the effective population size of WCA (as mirrored by the bottleneck signature detected in one of our analyses) occurring c. 4,300-7,440 ya, a time of climatic stability favourable to rainforest cover and that predated a major period of savannah expansion and forest fragmentation in western central Africa (starting c. 4000 ya).Rather than a paleoclimatic influence, we hypothesize that the arising of the Ceramic Late Stone Age culture in the region from c. 7000 BP onward -characterized by bifacial macrolithic and polished stone tools 40 -could have constituted the technological change (from "Stone to Metal Age"; 41 ) that made hunting more efficient and responsible for the decline of WBP.Such decline, although drastic (loss of c. 74.10-82.5% of the effective population size), places the effective populations size of WBP from WCA (2161 ± 219) within the conservative threshold range of minimum viable population size (500-5000).Because viable population size levels are paramount to maintaining acceptable levels of genetic diversity, it will be important to create or maintain forest corridors among WBP populations and mitigate deforestation and illegal hunting activities in order to preserve the genetic diversity levels of pangolin populations from the study region.

An endemic, trans-national trade of white-bellied pangolins from western central Africa, with Yaoundé as a main domestic hub
We showed that the 20 microsatellite loci specifically developed for WBP 21 provided the necessary power to confidently distinguish among all the WCA and Gab individuals.Only five microsatellite loci were needed to reach the conservative value of probability of identity < 0.01 42 , compared to the seven loci required for distinguishing pangolins from the Dahomey Gap 22 .This has important implications for the implementation of a convenient genetic tracing tool applied to the pangolin trade in western central Africa.Indeed, even a subset of our markers would be capable of estimating the exact number of individuals from scale seizures and trace the parallel trafficking networks of scale and body trade where same individuals can be dispatched 8 .Microsatellites genotyping also allowed identifying eleven cases where a same pangolin had been sampled twice, illustrating the difficulties of relying on third parties such as trained local assistants who are part of the bushmeat chain 14 , but at the same time showing the usefulness of our markers to counter such sampling bias.
Although mtDNA was not resolutive enough to trace the trade at the local scale, our results showed that the trade of WBP was endemic to western central Africa, involving both WCA and Gab lineages but excluding the other lineages described from central and West Africa.This is in line with the trade endemism also found for the Dahomey Gap lineage 22 , and suggests that large, urban bushmeat markets (here, Douala and Yaoundé for the most sampled) are not necessarily hubs for long-range, international trade of WBP.Rather, they represent domestic hubs that source large volumes of WBP at the sub-regional level (i.e.country-wide and neighboring countries), but are not directly involved in global trade networks such as those identified from harbor or airport seizures in Africa and Asia 9 .
The private allele frequency approach proposed by Zanvo et al. 22 allowed us to infer geographic origin hypotheses for 47 WBP sold from three urban bushmeat markets in Cameroon and seized at a European airport, representing c. 13% of the traceable individuals from WCA.Despite this relative success and the fact that our sample set covered quite exhaustively the WCA range, such inferences should be considered with caution as private allele frequency estimates will depend on the genetic representativeness of the sampled populations and the dispersal dynamics among populations 43,44 .Although we carefully controlled for the sampling effort (number of individuals) bias within populations using rarefaction, a denser geographic sampling of reference populations might dilute the tracing signal that we obtained in this study.In order to validate the hypotheses of trading sources that we generated, more powerful -genomic-approaches will be needed 23,26 .
Our results suggest that the Yaoundé market is the main national hub for the trade of bushmeat in Cameroon 27 , recruiting among six of the seven geographic sources that were traced as part of this study.Sourcing occurs in the western, southern and eastern part of the country, but also in Equatorial Guinea, representing a distance network www.nature.com/scientificreports/ between c. 120 and 370 km.Such inferences contrast with the more localized sources of the international trade that were recently inferred for the species 23 .
One geographic sourcing North of Douala and exclusive to the Yaoundé market, Manengolé, may provide evidence for transit of pangolins through Douala to Yaoundé, urban bushmeat markets in this case playing both a role of drop-off point and transit zone 27 .The two bushmeat markets from Douala showed contrasting spheres of influence, with Dacat being fed by a minimum of five geographic sources contra solely two sources for the Central market.Having said that, both showed a wide range of sourcing, respectively extending c. 600 and 500 km from Douala, beyond the maximal range delineated for Yaoundé.The two markets shared a unique supply source through the Yabassi population, situated c. 90 km north of Douala, in line with the recognized role of the unprotected Ebo forest as a feeding source of the Douala bushmeat market places 14 .
Eventually, we could trace the presence of WBP from the forests of southern Equatorial Guinea in Douala and Yaoundé markets, and Brussels airport (with flights originating from Cameroon).This result suggests that trans-national trade across neighbouring countries occur in the study region, potentially following the main road network connecting Equatorial Guinea to Cameroon, and that Cameroon constitutes a major hub of this sub-regional trade through the large feeding network of its urban bushmeat markets ( 8 ; https:// camer oonvo ice.com/ news/ 2021/ 07/ 22/ camer oun-trafi quants-de-pango lin-en-afriq ue-centr ale-demas ques/; this study).On the other hand, the likely occurrence of the Gab lineage and admixture between WCA and Gab or individuals from northern Gabon (Medouneu) in Cameroon blurred our tracing ability for a potential WBP trade occurring between Gabon and Cameroon, as previously mentioned in the literature 46 and recent media covers (https:// www.gabon media time.com/ 216-point es-divoi re-saisi es-camer oun-prove nant-gabon-react ions-de-lanpn-gouve rneme nt-atten dues/).
The number of haplotypes and number of private alleles were strikingly high in the Yaoundé market, supporting the hypothesis that Yaoundé is one of the main pangolin trade hubs of the study region.This was confirmed by the 12 different geographic sources mentioned by bushmeat vendors when asked about the origin of the pangolins sold in Yaoundé (N = 59; data not shown), including Sangmelima (and its surroundings), Mekin, Djoum, Akonolinga, Makenene, Bafia, Ma'an, Campo Ayos, Lomié, Yokadouma, Foumbot, Mbouda, Magba, Eseka, Mboumnyebel and Bankim.Such wide geographic network, obviously minimized by the genetic data set at hands, suggests a feeding network with a maximum sourcing distance of c. 600 km from the city, more in line with the large market's network range expected from the literature (see above).The great number of private alleles detected for the Yaoundé market (42 vs. 0 to 8 in reference populations and Douala markets, and 15 in Gabon) could suggest that, despite our considerable sampling effort, major reference populations remain unsampled.However, the relatively low private allele frequency calculated from Yaoundé supports the view that such high numbers of private alleles are an effect of sampling bias (> 250 samples).We believe that additional sampling of reference populations will have to be conducted across western central Africa to reach higher accuracy in the tracing of the WBP trade, although the theoretical numbers of sampled individuals for obtaining confident estimates of allele frequencies among populations for forensic purpose might be hardly achievable (see 48 ).

Conclusion
Our study provides an unprecedented assessment of the genetic status of WBP from western central Africa.Combining mtDNA sequencing and microsatellites genotyping of several hundreds of pangolins, we identified range overlap between the WCA and Gab mitochondrial lineages with limited introgression in southern Cameroon, so the two lineages can still be considered separate evolutionary units.Further sampling and sequencing efforts -particularly through high-throughput sequencing approaches-will have to be conducted across the study region, notably in Gabon where the number of genetic samples was low, in order to assess more accurately the admixture pattern between WCA and Gab and geographic structuring within WCA.Ideally, sampling will have to be oriented towards forest sites in order to minimize the potential bias of inaccurate sourcing radius of the local markets herein used as reference populations (notably relative to a potentially artifactual Wahlund effect).Such refined data will be crucial to establish a genetic-based management unit delineation strategy for guiding the conservation of WBP in the study region.As a corollary, new knowledge on the dispersal and breeding behavior of WBP, which remains literally unknown, will have to be produced urgently.
Our set of markers provides a unique opportunity to trace the WBP trade at the individual and population levels across western central Africa, and should be considered a valuable resource for future forensic applications at reasonable costs.For the geographic tracing of WBP from western central Africa, and because single locus (mtDNA)-typing has recently been applied to trace the pangolin trade 9 , we recommend a multi-locus approach given the admixture pattern observed between WCA and Gab.Because introgression may occur between other neighboring lineages of WBP 49 , it will be important in the future to rely on a global reference database where all the WBP lineages will have been genotyped.
According to our results and previous investigations, Yaoundé is a main hub for the pangolin trade in the sub-region, where WBP are openly sold on the stalls despite being nationally protected 14 .The Yaoundé market is driven by a large urban demand for bushmeat consumption, and as such should receive urgent attention from researchers, conservationists and anti-poaching services as a crucial point of entry to mitigate the WBP trade in Cameroon and, more globally, western central Africa.A similar conclusion applies to European points of entry of WBP carcasses, where CITES regulations apply but are in practice not efficiently implemented because of a global lack of prioritization, means and training at European borders 14 .As a general conclusion, better sensitization of range-state wildlife management services and foreign customs agencies to the regional and global trade of WBP should help prioritizing the conservation of the species.

Clustering analysis
We used cyt b to phylogenetically assign the collected samples to their respective mitochondrial lineages, according to Gaubert et al. (2016).For that purpose, we combined our samples from Cameroon (N = 505), Equatorial Guinea (N = 15), Bioko Island (N = 1), Gabon (N = 18) and European seizures (N = 12) to the cyt b sequences available for WBP at the time of the study (N = 76).A clustering distance tree was obtained with MEGA-X v10.2.2 using Neighbor Joining and the Kimura 2-parameter model 51 .Node support was assessed through 500 bootstrap replicates.

Genetic diversity and structure
Based on the new cyt b sequences that we produced, we reassessed the levels of genetic variation in WCA and Gab, after removing sequences with missing data (remaining N = 563).We used DnaSP v6.12 to compute number of polymorphic sites (S), haplotype number (h), haplotype diversity (Hd) and nucleotide diversity (π) (Supplementary Table S2 online).We mapped the geographic distribution of WCA and Gab haplotypes using ArcGIS 10.1 (Esri France).We ran Network v10.2.0.0 to build a median-joining haplotype network with ε = 0, so to minimize the presence of alternative median networks.

Demographic history
Because of the low number of samples attributed to Gab (N = 11; see Results), we explored the demographic history of WCA only (N = 552).We performed mismatch analysis using Arlequin 3.5.2 to test for signatures of demographic and spatial expansions, by calculating the sum of squared deviations (SSD) between observed and expected distributions using 1,000 bootstrap replicates.
Past variation in effective population size was assessed through a Bayesian skyline plot analysis in BEAST 2.6.3.We used the Bayesian Model Test option (bModelTest;) for site modeling, a coalescent Bayesian skyline prior and a relaxed clock (log normal-distributed).We fixed the mean mutation rate (ucld mean) to 0.026 substitution/site/Myrs with a normal distribution ranging from 0.010 to 0.045, as extracted from the posterior distribution values of the calibrated tree of Gaubert et al. 20 (data not shown).We ran 500.10 6 MCMC iterations, with trees and model parameters sampled every 10,000 generations and a 10% burnin.Parameter convergence and Bayesian skyline plot reconstruction were, respectively, visualized and performed in Tracer 1.7.1.Effective population size (N e ) was estimated considering a generation time of 2 years in WBP 22 .

Genetic diversity
We used Geneious 9.0.5 52 and the Microsatellites plugin (https:// www.genei ous.com/ featu res/ micro satel litegenot yping/) for allele scoring and genotype extraction.Only WCA and Gabon individuals with at least 75% of genotyping success were considered for the analyses (N = 558).The 75% threshold coincides with ≥ 15 microsatellite markers, which triples the minimum number of loci needed to discriminate among WCA individuals 21 53 , assuming population at equilibrium, on the Foumbot population (N = 31).Deviation from Hardy-Weinberg Equilibrum (HWE) was calculated in GenAlEx 6.5 54 .Linkage Deseliquilibrum (LD) was assessed using FSTAT 2.9.4 55 after 1000 randomizations (P-value ≥ 0.00005).The Bonferroni correction was applied to null hypothesis testing in those three analyses.
Genetic diversity was assessed through the number of alleles (Na), observed (Ho), expected (He) and unbiased expected (uHe) heterozygosity (in GenAlEx), allelic richness (A R ) and deficit of heterozygotes (F IS ) (in FSTAT), as mean values per population (see below).

Genetic structure
Global genetic variance among all genotyped individuals was assessed through pairwise population matrix unbiased genetic distances, and was visualized in GenAlEx using Principal Coordinates Analysis (PCoA).
For downstream population-based analyses, we defined the following 11 populations based on geographic proximity and biogeographic continuum among forests: Bayomen (Bayomen and Bangong  S1 online) were excluded because not proximate enough to any site to be merged and reach a sufficient number of samples.
We used STRU CTU RE 2.3.4 56 to conduct a clustering analysis, following a two-step dataset refinement procedure, where (i) all the individuals from the study area were included (WCA + Gab; N = 558), and (ii) only the WCA individuals were included (WCA; N = 181), discarding large urban markets and cyto-nuclear hybrids or admixed individuals (who with at least 20% of shared ancestry).The first scheme allowed screening the global genetic structure of WBP at the largest study scale.The second scheme aimed at assessing the genetic structure among WCA reference populations.The first and second schemes were run with the LOCPRIOR model, notably to manage the sampling bias among lineages (first scheme) and populations (second scheme; 57 ).We performed for each scheme 20 independent runs for K = 2 to 10 with 10 5 Markov chain Monte Carlo (MCMC) iterations (burnin = 10 4 ), assuming admixture and uncorrelated (first scheme) or correlated (second scheme) allele frequencies.We used STRU CTU RE HARVESTER 0.6.94 to estimate the most likely number of populations (K) using the ɅK method 58 .Graphic outputs from independent runs were summarized with CLUMPP and CLUMPAK server (http:// clump ak.tau.ac.il/) was used to generate the bar plots of population assignment.
We also inferred the most likely number of geographic populations using WCA georeferenced individuals, through the Geneland package 59 in R 4.0.5 (R Team Development Core 2021).We followed Coulon et al. 60 by first allowing K to vary from 1 to 20 and launched five runs of 5.10 5 MCMC iterations (thinning = 500; burnin = 500) under the frequency-correlated model and 1 km of uncertainty for spatial coordinates.Second, we fixed the number of estimated populations to K = 10-13 after the results of the first analysis [highest mean posterior densities for: K = 10 (0.21), K = 11 (0.30), K = 12 (0.26) and K = 13 (0.15)], and performed 20 independent runs with the same parameters.In the final analysis, the stability of individual assignment to populations was assessed among the best five runs (i.e. with the highest posterior probabilities).We used 500 × 500 pixels to map the posterior probabilities of population membership.
Pairwise differentiation (F ST ) among the 11 populations as defined above was computed in Arlequin 3.5 61 .Significant differentiation values (P ≤ 0.05) were estimated using 10 5 MCMC iteration chains and 10,000 dememorization steps.
We tested isolation-by-distance (IBD) within WCA individuals and among WCA populations through the R package pegas 62 using a Mantel test.The significance of the correlation (r) between individual-based genetic (Edward's) and geographic (Euclidean) distances was estimated through 10,000 permutations.

Tracing the pangolin trade
The discriminative power of our 20 microsatellite markers among WBP individuals was evaluated after three different approaches: (i) counting the number of identical genotypes among samples using the Multilocus tagging option in GenAlEx (suboption Matches), (ii) computing the probability of encountering the same genotype more than once by chance using the R package poppr (method = single; 63 ), and (iii) calculating the indices of unbiased probability of identity and probability of identity among siblings (uPI and PIsibs) using Gimlet 1.3.3 64.
Because we could not find a deep level of genetic structuring within WCA, we followed the protocol based on private allele frequency rarefaction described by Zanvo et al. 22 to trace the trade of WBP in our study area.We used the generalized rarefaction approach implemented in ADZE 65 to compute private allele frequencies (paf) among various combinations of populations from the 10 source populations as defined above (removing Gab).Because the original scheme of 10 populations yielded the greatest inferred values of paf (data not shown), we graphed the paf rarefaction curves (sample size rarefied from 2 to 5, 5 being the lowest number of individuals among populations) for each locus in these 10 populations.Only loci for which paf were > 40% and reached a plateau or showed an increasing trend at N = 5 for a given population were considered as potentially useful for tracing the origin of WBP found in the urban markets.We then cross-filtered these results with the private alleles observed in the actual populations (using GenAlEx), and eventually retained the loci that showed both high potential for tracing (ADZE output) and observed private alleles (GenAlEx output) for the given populations.www.nature.com/scientificreports/As a final step, we manually screened the private alleles present in WBP sampled from the Cameroonian large urban markets (Douala and Yaoundé) and international seizures to trace their geographical origins.

Demographic history
We tested for the signature of demographic bottlenecks in WCA (excluding hybrid and admixed individuals) using both the Single Mutation Model (SMM) and the Two Phase Model (TPM) in BOTTLENECK 1.2.02 66 .We applied the Wilcoxon sign-rank test to estimate heterozygote excess/deficit using 10,000 replications.
We assessed temporal changes in the effective population size of WCA through an approximate-likelihood approach implemented in the R package VarEff 67 .We followed Zanvo et al. 22 by fixing a generation time of 2 yrs and a microsatellites mutation rate of 5.10 -4 per generations 68 .We ran the analysis with three mutation models (single, geometric and two phases), using 10,000 MCMC batches (length = 1, thinned every 100 batches; JMAX = 3).The first 10,000 batches were discarded as burnin.

Figure 2 .
Figure 2. Bayesian skyline plots showing the demographic history of white-bellied pangolins from the Western Central Africa lineage.Y-axis = Effective population size (N e ); X-axis = Time in years before present.

Figure 3 .
Figure 3. Assignment plots among individuals of white-bellied pangolins from western central Africa as obtained with STRU CTU RE (left) and most likely number of population clusters (K = 2) following the ɅK method (right).Each individual is represented by a vertical bar.A-including all the samples from urban markets, seizures and reference populations (N = 558), sorted according to their mtDNA-based lineage assignment (WCA vs. Gab); asterisks indicate cyto-nuclear hybrids; X is an admixed individual (> 20% of genomic ancestry shared with Gab).B-only retaining the samples from 10 reference populations (N = 181); M = cluster from Medouneu (northern Gabon).

Figure 4 .
Figure 4. Temporal change in effective population size (Ne) of white-bellied pangolins from Western Central Africa, as estimated using VarEff under three different mutation models.Harmonic means (black line) and kernel density (color scale) of Ne posterior distributions are given in years BP.

No. private alleles (without urban markets) No. private alleles (with urban markets) Private allele frequency (without urban markets) Private allele frequency (with urban markets)
.Detection of null alleles in the 20 loci was performed with Microcheker 2.2.3